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Abstract: The velocity distribution function of dark matter particles is expected to show 
significant departures from a Maxwell-Boltzmann distribution. This can have profound ef- 
fects on the predicted dark matter - nucleon scattering rates in direct detection experiments, 
especially for dark matter models in which the scattering is sensitive to the high velocity 
tail of the distribution, such as inelastic dark matter (iDM) or light (few GeV) dark matter 
(LDM), and for experiments that require high energy recoil events, such as many direc- 
tionally sensitive experiments. Here we determine the velocity distribution functions from 
two of the highest resolution numerical simulations of Galactic dark matter structure (Via 
Lactea II and GHALO), and study the effects for these scenarios. For directional detection, 
we find that the observed departures from Maxwell-Boltzmann increase the contrast of the 
signal and change the typical direction of incoming DM particles. For iDM, the expected 
signals at direct detection experiments are changed dramatically: the annual modulation 
can be enhanced by more than a factor two, and the relative rates of DAM A compared 
to CDMS can change by an order of magnitude, while those compared to CRESST can 
change by a factor of two. The spectrum of the signal can also change dramatically, with 
many features arising due to substructure. For LDM the spectral effects are smaller, but 
changes do arise that improve the compatibility with existing experiments. We find that the 
phase of the modulation can depend upon energy, which would help discriminate against 
background should it be found. 
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1. Introduction 

Direct detection experiments aim to detect the low energy nuclear recoil from rare scatter- 
ing events between dark matter (hereafter DM) particles, assumed to be weakly interacting 
massive particles (WIMPs), and target nuclei. The event rate and its energy spectrum de- 
pend on the properties of the DM distribution at Earth's location, about 8.5 kpc from the 
Galactic Center [1]. Typical calculations of the scattering rate assume a "standard halo 
model" (SHM) consisting of a local DM density of 0.3 ± 0.1 GeV cm" 3 and a Maxwell- 
Boltzmann (MB) velocity distribution function with a (three-dimensional) dispersion of 
270 ± 70 km/s [2] truncated at an escape speed of ~ 550 km/s. Recent numerical simula- 
tions of the formation of Galactic-scale DM halos have reached the necessary resolution to 
directly test these assumptions. 

At the same time, absent a positive signal, a set of uniform halo assumptions allows 
a simple means to compare different experiments. However, in light of the recent results 
from DAMA/LIBRA [3], confirming earlier results from DAMA/Nal [4], it is important to 
consider halo model uncertainties when discussing exclusion limits from experiments with 
different targets, or energy ranges. This is particularly important because proposals such as 
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light dark matter (LDM) [5, 6] and inelastic dark matter (iDM) [7, 8], which aim to reconcile 
DAM A with null results from other experiments, sample the high velocity component of the 
WIMPs preferentially, and it is especially here that the Maxwell-Boltzmann distribution is 
expected to break down. 

In the hierarchical structure formation paradigm of standard cold dark matter cosmol- 
ogy the DM halo of a typical galaxy is built up through the merger of many individual 
gravitationally bound progenitor halos, which themselves were assembled in a hierarchical 
fashion. High resolution cosmological simulations, such as Via Lactea [9, 10], GHALO [11], 
and Aquarius [12], have shown that this merging process is in general incomplete, with 
the dense cores of many of the merging halos surviving as subhalos orbiting within their 
respective host halos. During pericenter passages tidal forces can strip off a large fraction 
of a subhalo's material, but the resulting cold tidal streams can readily be identified as 
velocity space substructure. The resulting DM halos are not perfectly phase mixed, and 
the assumption of a smooth halo is in general not a good one, neither in configuration 
space nor in velocity space [13]. 

At 8.5 kpc the Sun is located quite close to the Galactic Center, at least when compared 
to the overall extent of the Milky Way's DM halo (r v ; r ~ 200 — 300 kpc). This central 
region is notoriously difficult for cosmological numerical simulations to resolve, as very 
high particle numbers and very short time steps are required to avoid the so-called "over- 
merging problem" [14], which has until recently resulted in an artifically smooth central 
halo devoid of any substructure. With the advent of O(10 9 ) particle simulations at the 
Galactic scale this problem finally seems to have been overcome, with hundreds of subhalos 
identified at < 20 kpc. Nevertheless the local phase space structure is far from completely 
resolved, and likely never will be through direct numerical simulation. Any estimation of 
the importance of local density or velocity substructure based on cosmological simulation 
must thus rely on extrapolations over many orders of magnitude below its resolution limit. 

Local density variations due to the dumpiness of the DM halo are unlikely to signif- 
icantly affect the direct detection scattering rate. Based on the Aquarius Project suite 
of numerical simulations, Vogelsberger et al. (2008) [15] report that at more than 99.9% 
confidence the DM density at the Sun's location differs by less than 15% from the average 
over a constant density ellipsoidal shell. Extrapolating from their numerical convergence 
study they estimate a probability of 10~ 4 of the Sun residing in a bound subhalo of any 
mass. Analytical work by Kamionkowski & Koushiappas (2008) [16] predicts a positively 
skewed density distribution with local densities as low as one tenth the mean value, but 
probably not much less than half. 

The situation for velocity substructure is less clear. It is well established that nu- 
merically simulated dark matter halos exhibit significant velocity anisotropy and global 
departures from a Maxwell-Boltzmann distribution [17, 18, 19, 20]. The implications for di- 
rect detection experiments of these global departures from the standard Maxwellian model 
have previously been investigated in the context of a standard WIMP model [21, 22, 23, 24] 
and for inelastic dark matter [25, 26], and they were found to result in appreciable differ- 
ences (factor of a few) in the total event rates and the annual modulation signal. Most 
recently Vogelsberger et al. (2008) reported significant structure ("wiggles") in the ve- 
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locity distribution function measured in their high-resolution Aquarius simulations, which 
they attributed to events in the halo's mass assembly history. Their analysis concluded 
that velocity substructure due to bound subhalos or unbound tidal streams, however, does 
not influence the detector signals, since it makes up a highly sub-dominant mass fraction 
locally. 

The aim of this paper is take a closer look at this velocity space substructure and to 
examine its impact on the direct detection signal for models that are particularly sensitive 
to the high velocity tail, such as LDM or iDM. In contrast to Vogelsberger et al. (2008), we 
find that both global and local departures from the best-fit Maxwell-Boltzmann distribution 
can significantly affect the total event rate, the annual modulation, and the recoil energy 
spectrum. Parameter exclusion limits derived using a standard MB halo model are likely 
to be overly restrictive. 

This paper is organized as follows: in Section 2 we present velocity distribution func- 
tions derived from the high-resolution numerical simulations Via Lactea II and GHALO. 
In Section 3 we look at the implications of high velocity substructure for direct detection 
experiments with directional sensitivity. In Section 4 we consider iDM and LDM models 
and show how the observed local and global departures from the MB model affect scatter- 
ing event rates and recoil spectra at several ongoing direct detection experiments, and how 
this modifies parameter exclusion limits. A summary and discussion of our results can be 
found in Section 5. 

2. Results from Numerical Simulations 

The nuclear scattering event rate depends on the size of the detector, the type of target 
material, the scattering cross section, and the number density and velocity distribution of 
the impinging DM particles. We defer calculations of the expected event rate for various 
experimental setups and types of DM models to section 4, and focus in this section on the 
particle velocity distributions, which we determine directly from numerical simulations. 

2.1 The Via Lactea and GHALO simulations 

Our analysis is based on two of the currently highest resolution numerical simulations of 
Galactic DM structure: Via Lactea II (VL2) [10] and GHALO [11]. Both are cosmological 
cold DM N-body simulations that follow the hierarchical growth and evolution of a Milky- 
Way-scale halo and its substructure from initial conditions in the linear regime {z = 104 
for VL2, z = 58 for GHALO) down to the present epoch. For details about the setup of the 
simulations we refer the reader to the above references. The VL2 host halo is resolved with 
~ 400 million particles of mass m p = 4, 100 M Q within its virial radius 1 of r v - a = 309 kpc 
and has a mass of Mh a i = 1.7 x 10 12 M Q and peak circular velocity V max = 201.3 km/s. 
The GHALO host is somewhat less massive, Mh a i Q = 1.1 x 10 12 M Q and V max = 152.7 
km/s, but even more highly resolved, with 1.1 billion particles of mass m p = 1,000Mq 
within its r v ; r = 267 kpc. For reference we show the circular velocity of the two halos in 

Vvir is defined as the radius enclosing a density of A v i r = 389 times the background density. 
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Figure 1: Circular velocity profiles of the VL2, GHALO, and GHALO s host halos. 

Fig.l. In order to facilitate a more direct comparison between the two halos, we have also 
scaled GHALO to match VL2's V mSuX by multiplying the simulation's length and velocity 
units by a factor / = V r max (VL2)/V max (GHALO) = 1.32, and the mass unit by / 3 . We 
refer to this model as GHALO s . The circular velocity of these three halos at 8.5 kpc is 
158.1, 121.7, and 148.9 km/s for VL2, GHALO, and GHALO s , respectively. 

2.2 Velocity Modulus Distributions 

The DM-nucleon scattering event rate is directly proportional to 



where f(v) is the DM velocity distribution function in the Earth's rest frame and v m i n (En) 
is the minimum velocity that can result in a scattering with a given nuclear recoil energy Er. 
For a target with nuclear mass and a WIMP/nucleon reduced mass fi = m]^m x /(mN + 
m x ), v m i n (E R ) is given by 



The 5 refers to the possible mass splitting between the incoming and outgoing DM particle, 
which would be for standard and light DM and 0(100 keV) for inelastic DM. 

We determine f{v) in the halo rest frame directly from the particle velocities in our 
numerical simulations, and in the Earth's rest frame by first applying a Galilean velocity 
boost by v©(t). The Earth's velocity with respect to the Galactic center is the sum of the 
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local standard of rest (LSR) circular velocity around the Galactic center, the Sun's peculiar 
motion with respect to the LSR, and the Earth's orbital velocity with respect to the Sun, 

v®(t) = v LSR + v pec + v rut(t)- (2-3) 

We follow the prescription given in Chang et al. (2008) [27] and set vlsr = (0, 220, 0) 
km/s, iTpec = (10.00,5.23,7.17) km/s [28], and v OT bit(t) as specified in reference [29]. The 
velocities are given in the conventional (U, V, W) coordinate system where U refers to 
motion radially inwards towards the Galactic center, V in the direction of Galactic rotation, 
and W vertically upwards out of the plane of the disk. We associate these three velocity 
coordinates with the (v r ,vg, v^) coordinates of the simulation particles. 

The Earth's orbital motion around the Sun results in the well-known annual modu- 
lation of the scattering rate, which the DAM A collaboration claims to have detected at 
very high statistical significance [3]. For the SHM the peak of this modulation occurs 
around June 2 nd , when the Earth's relative motion with respect to the Galactic DM halo 
is maximized. 

We have measured the DM velocity distribution from all particles in a 1 kpc wide 
spherical shell (8 kpc < r < 9 kpc), containing 2.1, 5.4, and 3.6 million particles in VL2, 
GHALO, and GHALO s , respectively. The large particle numbers in these measurements 
result in a very small statistical uncertainty, but fail to capture any local variations. To 
address this we have also determined f(v) from the particles in 100 randomly distributed 
sample spheres centered at 8.5 kpc. These sample spheres have radii of 1.5 kpc for VL2 
and 1 kpc for GHALO and GHALO s , and contain a median of 31,281, 21,740, and 14,437 
particles in the three simulations. 2 

The resulting distributions, both in the halo rest frame and translated into Earth's rest 
frame, are shown in Fig. 2. The shell averaged distribution is plotted with a solid line, while 
the light and dark green shaded regions indicate the 68% scatter around the median and the 
absolute minimum and maximum values of the distribution over the 100 sample spheres. 
For comparison we have also overplotted the best-fitting Maxwell-Boltzmann (hereafter 
MB) distributions, with ID velocity dispersion of <7id = 130, 100, and 130 km/s. These 
clearly underpredict both the low and high velocity tails of the actual distribution. This 
is not a new result and has previously been found in cosmological numerical simulations 
[17, 18, 19, 15]. Actually there is no reason to assume that a self-gravitating, dissipationless 
system would have a locally Maxwellian velocity distribution, and in fact it has been shown 
that self-consistent, stable models of cuspy DM structures require just such non-Gaussianity 
[30, 31]. 

In addition to its overall non-Maxwellian nature, we notice several broad bumps present 
in both the shell averaged and, at very similar speeds, in the sub-sample f{v). Similar fea- 
tures were reported by Vogelsberger et al. (2008) [15] for the host halos of their completely 
independent Aquarius simulations, and thus appear to be robust predictions of hierarchi- 
cally formed collisionless objects. Vogelsberger et al. also showed that the broad bumps 

2 Tables of <?(u m i n ) determined from the spherical shell and the 100 sample spheres, and trac- 
ing the annual modulation over 12 evenly spaced output times, are available for download at 
http : //astro .berkeley . edu/~mqk/dmdd/. 
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Figure 2: Velocity distribution functions: the left panels are in the host halo's restframe, the 
right panels in the restframe of the Earth on June 2 nd , the peak of the Earth's velocity relative 
to Galactic DM halo. The solid red line is the distribution for all particles in a 1 kpc wide shell 
centered at 8.5 kpc, the light and dark green shaded regions denote the 68% scatter around the 
median and the minimum and maximum values over the 100 sample spheres, and the dotted line 
represents the best-fitting Maxwell-Boltzmann distribution. 



are independent of location and persistent in time and hence reflect the detailed assembly 
history of the host halo, rather than individual streams or subhalos. The extrema of the 
sub-sample distributions, however, exhibit numerous distinctive narrow spikes at certain 
velocities, and these are due to just such discrete structures. Note that although only 
a small fraction of sample spheres exhibits such spikes, they are clearly present in some 
spheres in all three simulations. The Galilean transform into the Earth's rest frame washes 
out most of the broad bumps, but the spikes remain visible, especially in the high veloc- 
ity tails, where they can profoundly affect the scattering rates for inelastic and light DM 
models (see Section 4). 
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In order to assess the dependence of these features on the sample sphere size, we also 
considered for the VL2 simulation sphere radii of 1 and 2 kpc, containing a median of 
9,200 and 74,398 particles, respectively. The coarse features of the distributions persist, 
but of course the prevalence of the spikes increases with sample sphere size, as more of the 
substructure is probed. It is difficult to assess with our simulations the true likelihood of 
significant local velocity substructure, as it depends on the abundance and physical extent 
of subhalos and tidal streams many orders of magnitude below the length scales that we can 
accurately resolve. Higher resolution numerical simulations, as well as analytical models 
[32], perhaps in conjunction with simulations [33], will be necessary to settle this question. 

2.3 The effects of neglected baryonic physics 

The values of <tid we report here may appear surprisingly low to a reader familiar with the 
standard isothermal MB halo assumption of (v 2 ) = 3o"^ D = 3/2 Vq, where vq, the peak of 
the MB distribution, i.e. the most probable speed, is assumed to be equal to the rotation 
velocity of the Sun around the galaxy, vq ~ 220 km/s. In this standard model ctid would 
be 156 km/s, considerably higher than our values of 130 km/s and 100 km/s, respectively. 
In fact, the local circular velocity v c and velocity dispersion a are only indirectly related 
and not necessarily equal. While v c is set by the local radial gradient of the potential, a 
depends on the shape of the potential at exterior radii. For a non-rotating spherical system 
the relation between v c {r) and the radial velocity dispersion <r r (r) is given by 

2 2 / dlnp dlna 2 \ 

" c = -^idIn7 + 7Iny + 2/3 J' {2A) 

where = 7(7") is the logarithmic slope of the density profile and (3{r) = 1 — a^/a 2 is 
the velocity anisotropy. In a singular isothermal sphere (7 = —2, d J" n ^, r = 0, (3 = 0) we 
have v c = vq, but for an NFW profile v c /v ~ 0.88 at r = r s /2. In the VL2 and GHALO 
host halos the relation is v c /vq = 0.85 and 0.86, respectively. 

A central baryonic condensation in the form of a Galactic disk and bulge will deepen the 
central potential, raise the local circular velocity to ~ 220 km/s, and increase the velocity 
dispersion of DM particles at the Sun's location. Since the VL2 and GHALO simulations 
do not include baryonic physics, it is not surprising that the values of v c and vq at 8.5 kpc 
are lower than appropriate for our Milky Way galaxy. The main focus of our work here 
is to investigate the effects of global and local variations from the MB assumption, and 
therefore we compare our results to the best-fitting MB model (vq = 184 km/s for VL2 and 
GHALO s and 141 km/s for GHALO) instead of the standard MB halo model (v = 220 
km/s). In principle we could have scaled just the velocities up to give vo = 220 km/s, 
but this would remove many of the effects of substructure by pushing it above the escape 
velocity. Thus, we use VL2 and GHALO s simulations for the parameter exclusions plots 
in Section 4, but it should be recognized, that the velocity dispersion is somewhat lower 
than in conventional MB halo parameterizations. 

2.4 Local Escape Speed 

The escape speed from the Sun's location in the Galactic halo is another factor that can 
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Figure 3: Radial and tangential velocity distribution functions. The solid lines show the shell 
average, the shaded range the 68% scatter around the median, and the dotted curve shows the 
best-fit MB model. 

strongly affect scattering rates, especially for inelastic and light DM models. The RAVE 
survey's sample of high- velocity stars constrains the Galactic escape velocity to lie between 
498 and 608 km/s at 90% confidence, with a median likelihood of 544 km/s [34]. This is 
in good agreement with the highest halo rest frame speed of any particle in our 8.5 kpc 
spherical shells, namely 550 km/s in VL2 and 586 km/s in GHALO s . The lower V max of 
the GHALO host is reflected in a significantly lower escape speed, only 433 km/s. The 
corresponding maximum speeds in the Earth rest frame are 735-761, 773-802, and 634-660 
km/s, where the range refers to the modulation introduced by the Earth's orbit around 
the Sun. 

2.5 Radial and Tangential Distributions 

For comparison with past work [25, 24] we also present separately the distribution functions 
of the radial v r and tangential vt = ^Jvg + velocity components in Fig. 3. We fit these 
distributions with functions like the ones in Fairbairn & Schwetz (2009) [24], except that 
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we don't normalize the velocities by the square root of the gravitational potential at the 
particles' location and instead provide separate fits for each simulation. We also include 
an estimate of the variance due to local velocity substructure by separately fitting the 
distribution in each of the sample spheres. The fitting functions are 



f(Vr) 
f{vt) 



1 




Nr 


exp 


Vt 




N t 


exp 



(2.5) 
(2.6) 



and the parameters of these fits are listed in Table 2.5. The normalizations N r and Nt 
can readily be obtained numerically for a given set of parameters by ensuring that the 
distributions integrate to unity. 
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148.9 


Oi r ,t 


0.934 


0.941 


0.877 


0.985 


0.642 


0.657 


0.638 


0.674 


GHALO 


v T)t [km/s] 


167.9 


163.6 


156.4 


173.0 


103.1 


114.3 


93.21 


137.0 


a T ,t 


1.12 


1.11 


1.02 


1.20 


0.685 


0.719 


0.666 


0.819 


GHALO s 


v T)t [km/s] 


217.9 


213.8 


202.3 


226.6 


138.2 


162.2 


125.1 


183.1 


a r ,t 


1.11 


1.11 


1.01 


1.18 


0.687 


0.759 


0.664 


0.842 



Table 1: Radial and tangential velocity distribution fit parameters (see Eqs. 2.5 and 2.6). The 
columns labeled "shell" refer to the fit for all particles in the 1 kpc wide shell centered on 8.5 kpc, 
whereas the following three columns give the median as well as the 16 th and 84 th percentile of the 
distribution of fits over the 100 sample spheres (see text for more detail). 



2.6 Modulation Amplitude and Peak Day 

The annual modulation of the scattering rate <?(i> m in) (Eq. 2.1) grows with v m i n , since a 
reduction in the number of particles able to scatter makes the summer-to- winter difference 
relatively more important. At sufficiently high velocities the modulation amplitude can 
even reach unity, when during the winter there simply aren't any particles with velocities 
above u m i n . 

The expected modulation amplitude for a given experiment is then determined by the 
relation between the measured recoil energy Er and v m i n , see Eq. 2.2. There are important 
qualitative difference between standard DM models (5 = 0) and inelastic models (<5 ~ 100 
keV). Inelastic DM models require a higher v m { n for a given Er, resulting in a lower total 
scattering rate and a more pronounced modulation. Furthermore, while for standard DM 
Vmm grows with Er, in inelastic models v m i n typically falls with Er (for 5 > Ertun / 7-0 • 

In the left panels of Fig. 4 we show the u m i n dependence of the annual modulation 
amplitude, (max(g(v m in)) ~ min)/(max + min) as measured in our simulations. We focus 
on the high v m i n region that is relevant for inelastic and light dark matter models, as well 
as directional detection experiments. In VL2 and GHALO s the shell averaged modulation 
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Figure 4: The amplitude (left panels) and peak day (right) as a function of w m i n . The solid red 
line indicates the shell averaged quantities, and the dotted line the best-fit MB distribution. The 
light and green shaded regions cover the central 68% region around the median and the minimum 
and maximum values of the distribution over the 100 sample spheres. The thin black line shows 
the behaviour of one example sample sphere. 

amplitude (solid red line) rises from about 20% a/t t'min 

= 400 km/s to unity at ~ 750 
km/s. The GHALO amplitudes are shifted to lower i> m i n , growing from 40% at 400 km/s 
to 100% already at ~ 600 km/s. The strong high velocity tails of f(v) (Fig. 2) result in 
somewhat lower modulation amplitudes compared to the best-fit MB distributions (dotted 
line). At the very highest velocities f{v) drops below the Maxwellian distribution in VL2 
and GHALO, and this leads to the rise in amplitudes above the Maxwellian case for v m \ n > 
600 and 550 km/s, respectively. In GHALO s the distribution more closely follows the MB 
fit, and only barely rises above it at t> m i n > 670 km/s. Interestingly, all three simulations 
exhibit a pronounced dip in the modulation amplitude at close to the highest v m i n . These 
correspond to bumps in f(v) discussed in Section 2.2. 

As before, the light and dark shaded green regions in Fig. 4 cover the 68% region 
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around the median and the extrema of the distribution over the 100 sample spheres. Over 
most of the range of t> m i n , the typical modulation amplitude in a sample sphere differs from 
the spherical shell average by less than 10%. The extrema of the distribution, however, 
can differ much more significantly, especially at higher velocities. If the Sun happens to be 
passing through a fast moving subhalo or tidal stream, the typical velocity of impinging 
DM particles can greatly differ from the smooth halo expectation, leading to an increase 
(or decrease) in the overall scattering rate and modulation amplitude. For a conventional 
massive (> 10 GeV) DM particle, scattering events are dominated by the peak of the 
velocity distribution function and the effects from streams or subhalos are washed out and 
typically negligible [15]. However, whenever scattering events are dominated by particles 
on the high velocity tail of the distribution, either because of a velocity threshold (inelastic 
DM), a particularly low particle mass (light DM), or for detectors intrinsically sensitive 
to only high recoil energies (e.g. many directionally sensitive detectors), the presence of 
velocity substructure can have a profound impact on the event rates. 

Another potentially interesting signature of velocity space structure is the shift in 
the peak day of the annual modulation, as was explored in [22] . For an isotropic velocity 
distribution (e.g. MB) the peak day is independent of w m i n and occurs around the beginning 
of June. Any kind of departure from isotropy, however, would leave its signature as a 
change in the phase of the modulation. In the most extreme case of a very massive DM 
stream moving in exactly the same direction as the Sun, the phase of modulation could 
flip completely, with the maximum scattering amplitude occurring in the winter. 

In the right panel of Fig. 4 WG sllOW tllG V min 

dependence of the peak day as measured 
in our simulations. The spherical shell average is consistent with a phase shift of zero, 
except at the very highest velocities where a few discrete velocity structures dominate 
the shell average and introduce small velocity dependent phase shifts. The peak days 
in the sample spheres, however, often differ by ~ 20 days from the fiducial value. As an 
example we have plotted a curve for one of the sample spheres. The peak day determination 
appears to be quite noisy here, but since the particle numbers are still fairly large (iV(> 
v min ) = 9114,802,105 for v min = 400,600,670 km/s in the VL2 sample sphere shown), 
we believe that the variations in peak day arise from actual velocity structure rather than 
discreteness noise. The 68% scatter of the phase shifts over the 100 sample spheres remains 
remarkable constant at ±10 — 20 days throught most of the t> m in range. At higher velocities 
(vmin > 600 — 700 km/s) the typical phase shifts grow to ~ 50 days, which could be 
due to the increasing relative importance of individual streams, but may also be due to 
particle discreteness noise. Based on our measurements we would expect some amount 
of variation in the peak day as a function recoil energy. As it is hard to imagine any 
background contamination to exhibit such a phase shift, this raises the possibility of such a 
measurement confirming a DM origin of the DAMA modulation signal. 

We conclude this section by noting that our simulations reveal both global (shell av- 
eraged) and local (sample spheres) departures from the standard halo model in velocity 
space. These can have a significant effects on the overall scattering rate, as well as on 
the amplitude and peak day of the annual modulation thereof. In the next section we go 
on to explore how these effects translate into predicted detection rates for actual direct 



- 11 - 



detection experiments, and how they modify parameter exclusions plots based on existing 
null-detections. 



3. Velocity Substructure and Directional Detection 

Up to this point we have focused on the integrated features of the halo, i.e., f(v) or 
ff(^min)) but there are many more structures that appear that are not pronounced in these 
measures. In particular, the presence of clumps and streams, while contributing to the 
bumps and wiggles of these functions, can be more pronounced when the direction of 
the particles' motions, rather than just the amplitude, is considered. These features are 
especially important for directional WIMP detectors, such as DRIFT [35], NEWAGE [36] 
and DMTPC [37], where the presence of such structures could show up as a dramatic 
signal. Although we do not discuss it here, such directional experiments have been argued 
to be especially important for testing inelastic models [38, 39]. 

To quantify this, we begin by searching for "hotspots" on the sky. To do this, we make 
a map of the sky in HealPix, and consider the flux of WIMPs from each direction in the 
sky. We divide the sky into 192 equal regions of 215 square degrees, and determine pi, 
the fraction of particles above v m i n with a velocity vector pointing towards bin i. We take 
the same 100 sample spheres as before, but note that beyond determining sample sphere 
membership we do not consider the location of a particle: the assignment to a given sky 
pixel is based solely on the direction of its velocity vector. All of the structure in a given 
sample sphere is considered local, i.e. able to influence the signal at an Earth-bound direct 
detection experiment. 

As an example of the kind of effects high velocity substructure can produce, we show 
in Fig. 5 the speed distribution f(v) and skymaps for VL2 sample sphere #03. In the 
halo rest frame (top row, v m \ n = 400 km/s) a very pronounced feature is visible due to the 
presence of a subhalo moving with galacto-centric velocity modulus of ~ 440 km/s. This 
feature persists in the Earth rest frame (center row, v m ; n = 500 km/s), where the direction 
of the subhalo's motion is "hotter" than the "DM headwind" hotspot in the direction of 
Earth's motion. When translating to the Earth's rest frame, there is an additional degree 
of freedom arising from the unspecified plane of the Galactic disk. We associate radial 
motion towards (away from) the Galactic Center with the center (anti-center) of the map, 
but the hotspot arising from Earth's motion is free to be rotated around this axis. In this 
example we have chosen this rotation to maximize the angle between the subhalo hotspot 
and Earth's motion. In the bottom row we show for comparison the Earth rest frame map 
for the MB-case without any substructure. 

This is just one strong example, and we would like to understand what sorts of hotspots 
we might expect. For this purpose we define the hotspot ratio 

irp / \ max fo(-tw)} 

the ratio of the hottest pixel in the sphere sample skymap above v m i n to the hottest pixel 
in the corresponding MB-case, and the hotspot angle tp between these two pixels. We 



- 12 - 




Earth Rest Frame (Maxwell— Boltzmann) 




300 400 600 

v [km/s] 




Figure 5: f(v) and HealPix skymaps of the fraction of particles above v m j n coming from a given 
direction, for VL2 sample sphere #03 which contains a fast-moving subhalo. The top row is 
in the halo rest frame (u m i n = 400 km/s), and the middle translated into the Earth rest frame 
{vmin = 500 km/s). For comparison the bottom row shows the Maxwell-Boltzmann halo case 
without substructure. 



calculate HR and ip for all 100 sample spheres, and in each sphere for a full 2tt rotation (in 
one degree increments) of the direction of Earth's motion. We show in Fig. 6 for VL2 the 
distribution of HR as a function of v m [ n and the distribution of ip for the case without a 
velocity threshold and for v m [ n = 500 km/s. For small velocities the mean of HR is unity 
and the r.m.s. variation is only 10%. As u m ; n increases, HR grows: at v m \ n = 600 km/s, 
the mean HR is 1.3±0.35, and at v m \ n = 700 km/s it's 3.1 ± 1.4. The downturn at the very 
highest velocities is caused by running out of particles in the sample spheres. There are 
also marked changes in the direction of the hottest pixel. Even without a velocity threshold 
{vram = km/s) in 38% of all cases the direction of the hottest spot on the sky is more 
than 10° removed from the direction of Earth's motion (i.e. the MB halo expectation). At 
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Figure 6: The distribution of HR(v m i„) as a function of w m i n (left) and of ip for v m i n — and 500 
km/s (right) for the VL2 simulation. 

^rnin = 500 km/s t/ie fruffi; of DM particles are more likely to be coming from a local hotspot 
with V> > 10° than from the direction of Earth's motion. In this case there is only an 18% 
chance of having tfj < 10°. Not all of this structure is due to individual subhalos; some of 
it may arise from tidal streams, and some from the anisotropic velocity distribution of the 
host halo particles. 

4. Implications for DAMA 

In light of the DAMA/LIBRA results [3], two scenarios which have garnered a great deal 
of attention of late are light dark matter [5, 6], and inelastic dark matter [7, 8]. These 
scenarios were proposed some time ago to address the conflicts between DAMA [4] and 
CDMS [40] as well as other experiments at the time. 

Recent studies [27, 25, 41, 42] have shown that iDM remains a viable explanation of the 
DAMA data, consistent with recent results from CDMS [43], XENON10 [44], KIMS [45], 
ZEPLIN II [46], ZEPLIN III [47] and CRESST [48]. Such models have some tension with 
CRESST (see the discussion in [27, 49]), which observed 7 events in the tungsten band, 
while approximately 13 would be expected [27]. While lighter masses have no tension with 
CDMS, higher mass (> 250 GeV) WIMPs do. 

Light dark matter no longer works in its original incarnation [5, 6], but instead relies 
upon "channeling," [50, 51], a difficult to quantify effect whereby some fraction of nuclear 
recoils have essentially all of the energy converted into scintillation light, rather than just 
a fraction. Even including channeling, such scenarios seem to have strong tension with 
the data [52, 24, 53], both in the spectrum of the modulation, and constraints from the 
unmodulated event rate at low (1-2 keVee) energies 3 . If one disregards the lowest (2-2.5 

3 keVee stands for "electron equivalent keV" , a unit of energy used for scintillation light, which is produced 
by interactions of the recoiling nucleus with electrons. It is related to the full nuclear recoil energy (in keVr) 
through the quenching factor q: £(scintillation)/keVee = q i?(recoil)/keVr. q is a material-dependent 
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keVee) bin from DAMA, however, the fit improves [52, 54], but there is still tension for the 
lightest particles from the presence of modulation above 4 keVee and an overprediction of an 
unmodulated signal at lower energies. It is important to note that the uncertainties in L e Q , 
the scintillation efficiency of liquid Xenon (see [55] for a discussion), are not completely 
included in these analyses, and the most recent measurements of L c g [55] suggest that 
the low energy analysis threshold may be somewhat higher than the 4.5keVr used with 
-^eff = 0.19, weakening the limits. On the other hand, the most recent analysis from 
XENON10 [56] using a more advanced rejection of double-scatter events shows no events 
at all in the signal region all the way down to the S2 threshold, strengthening the limits. 
In light of these uncertainties, and in view of the importance of the result, we believe it is 
best to maintain an open mind about the viability of the light WIMP explanation. 

For our purposes, the crucial feature is that these scenarios sample the high (v ~ 
600 km/s) component of the velocity distribution, and so are especially sensitive to the 
departures from a simple Maxwell-Boltzmann distribution. Consequently, it is important 
to investigate whether these changes affect the tensions described above. 

4.1 Inelastic Dark Matter 

The inelastic Dark Matter scenario [7] was proposed to explain the origin of the DAMA 
annual modulation signal. The scenario relies upon three basic elements: first, the presence 
of an excited state x* °f the dark matter \i with a splitting 5 = m x * — m x ~ m x v 2 
comparable to the kinetic energy of the WIMP. Second, the absence of, or at most a 
suppressed elastic scattering cross section off of nuclei, i.e., the process xiV — ► %iV should 
be small. Third, an allowed cross section for inelastic scatterings, i.e., %iV — ► x*iV, with a 
size set roughly by the weak scale. 

Although these properties may seem odd at first blush, they are in fact perfectly natural 
if the scattering occurs through a gauge interaction [7, 8] , where the splitting is between 
the two Majorana components of a massive pseudo-Dirac fermion, or between the real and 
imaginary components of a complex scalar. Simple models can be constructed where the 
mediating interaction is the Z-boson [7, 8, 57, 41]. Models with new vector interactions to 
explain the PAMELA positron excess also naturally provide models of iDM [58], and often 
where all scales arise naturally from radiative corrections [58, 59, 60, 61, 62]. Composite 
models provide a simple origin for the excited states as well [63, 64]. 

The principle change is a kinematical requirement on the scattering, 



where is the mass of the target nucleus and /i is the reduced mass of the WIMP/target 
nucleus system. For 5 ~ /i/3 2 /2, the consequences can be significant. This requirement 
has three principal effects: first, the kinematical constraint depends on the target nucleus 
mass, and is more stringent for lighter nuclei. If we are sampling dominantly the tail of the 
velocity distribution, the relative effect between heavy and light targets (e.g., iodine versus 

quantity that must be experimentally determined. 




(4.1) 



- 15 - 



germanium) can be significant. Second, again because the signal is sampling the tail of 
the velocity distribution, the modulated amplitude can be significantly higher than the few 
percent expected for a standard WIMP. In fact, in the cases where there are particles kine- 
matically scattering in the summer, but not the winter, the modulation can reach 100%. 
Third, the inelasticity suppresses or eliminates the event rate at lower energies. Because 
standard WIMPs have exponentially more events at lower energies, most experiments have 
focused on controlling background in this region and lowering the threshold. In partic- 
ular, the XENON10 and ZEPLINIII experiments have few events at low energies, but a 
significant background at higher energies. Consequently, their limits for standard WIMPs 
are quite strong, but significantly weaker for inelastic WIMPs. The combination of these 
three elements allows an explanation of the DAMA result that is consistent with all other 
current experiments [27, 25, 41]. 

4.1.1 Fitting the DAMA signal 

DAMA reports an annual modulation in the range of 2 — 6 keVee range. This can be 
interpreted as arising from either Na or I scattering events. In the former case (Na), it 
has been shown that the "light inelastic" region (an approximately 15 GeV WIMP with 
5 ~ 30 keV) can open up significant parameter space [52, 49]. In the latter case (I), we 
find the "heavy inelastic" region (a 100+ GeV WIMP with 5 > 100 keV) opens significant 
parameter space. Since the constraints are stronger on the heavy case, we focus on this 
region. 

When determining the precise values of parameters that might agree with DAMA, 
we must convert from keVee to keVr, which already introduces significant uncertainties, a 
point which has been recently discussed by [25]. The quenching factor for iodine has been 
found to have a range of different values, including q = 0.09 ±0.01 [65], q = 0.05 ±0.02 [66], 
q = 0.08±0.02 [67] and q = 0.08±0.01 [68] 4 . The first two measurements of the quenching 
factor are somewhat indirect, fitting event distributions at low energies. The last two are 
more direct. Of the last two, we should note that the first includes non-linearity in its 
stated uncertainty, while the second explicitly fixes to a signal electron energy, and thus 
this effect is not included. We adopt q = 0.08 ± 0.02 as the value for quenching, which 
corresponds to a range of 25 — 75 keVr for DAMA. 

4.1.2 Constraining the iDM interpretation of DAMA 

While iDM allows DAMA to exist consistently with the other experiments, the reach of 
other experiments is at or near the predicted DAMA level. In particular, as we see in 
Fig. 14, within the context of a MB halo, at high masses, CDMS excludes the entire range 
of the parameter space that fits DAMA. At the same time, the CRESST results (with 
a tungsten target) seem very tense with the data over the whole mass range. These two 
experiments provide the greatest present tension with the iDM interpretation of DAMA. An 

4 The value quoted by [68] is generally q = 0.086 ± 0.007, but this value averages over a wide range 
of energies. The two measured values for the DAMA energy range specifically are q = 0.08 ± 0.01 and 
q = 0.08 ±0.02. 
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Figure 7: The ratio of DAMA signal in a given simulation to that in a MB distribution. The 
solid red line is the spherical shell average, the dashed line the median of the distribution over the 
sample spheres. 



immediate question we can ask is then whether the velocity-space distribution of particles 
in the simulations makes the constraints more or less significant than a MB distribution. 

This is not an easy question to answer. Since there can be significant spatial variation 
in the velocity distribution, we would like to quantify this without making full exclusion 
plots for every data point. The three principle constraints on the iDM interpretation 
of DAMA come from a) Xe experiments (XENON10, ZEPLIN-II and ZEPLIN-III), b) 
CDMS (Germanium) and c) CRESST (Tungsten). While varying the halo model can have 
significant effects on each experiment, including rates and spectra, the variation of velocity 
structures in the different halos affects these limits differently, and it is difficult to quantify 
in aggregate what the effects are. 

Before proceeding into a detailed analysis of the DAMA signal, we can already study 
a preliminary question: how does the cross section needed to explain DAMA compare 
between a MB halo and a simulation? We consider the ratio of the modulation between 
MB and the simulations in Fig. 7. Overall, we see that the signal at DAMA is increased 
in the simulation, as much as a factor of a few. 

Of course, increasing the DAMA signal does not change the constraints if the sig- 
nal in any of the other experiments is changed, as well. We proceed by next consider- 
ing a ratio of ratios (RoR). Since Germanium (CDMS) has a higher velocity threshold 
than Iodine (DAMA), while Tungsten (CRESST) has a lower velocity threshold, a sim- 
ple test is to look at the variation of the signal at different experiments as a fraction of 
the DAMA modulation amplitude. We thus calculate for each of the simulation samples 
(spherical shell and 100 spheres), as well as for the best-fitting MB model, the ratios of 
the CDMS and CRESST signals, integrated from 10-100 keV, to the total DAMA mod- 
ulation in either of the high-q and low-g benchmark regions described above. We then 
divide the simulation ratio by the MB ratio: (CDMS / [DAMA] ) s ; m / ( CDMS/ [DAMA] ) mb 
and (CRESST / [DAMA]) s ; m / (CRESST / [DAMA] ) mb • If these RoR's are smaller than one, 
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Figure 8: Scatter plots of the ratios of ratios (RoR, see text for details) for the CRESST and 
CDMS experiments. RoR's less than one indicate that using the simulation's velocity distribution 
instead of the best-fitting MB model weakens the CRESST or CDMS limits relative to the DAMA 
(high-g) signal. 

it suggests that using the simulation's velocity distribution instead of the best-fitting MB 
distribution weakens the limits relative to DAMA, while values larger than one imply 
stronger limits. The effects on Xenon limits are harder to quantify, because of the added 
uncertainties in the conversion from keVee to keVr at those experiments, and where, pre- 
cisely, the backgrounds lie. Thus, we consider first only CDMS and CRESST limits. 

In Fig. 8 we show scatter plots of the CRESST RoR against the CDMS RoR. Large 
filled symbols indicate the spherical shell sample and small symbols are used for the sample 
spheres. Note that in some cases the CDMS RoR is zero, indicating that the simulation 
velocity distribution resulted in no CDMS signal at all. A few things are immediately 
obvious from this plot. First, the limits from CDMS can vary wildly between simulations, 
and even between different spheres within a single simulation. This simply represents how 
dramatically the velocity distribution can change at the highest velocities. Second, we see 
that the CRESST rate is much more weakly affected, with typically suppressions of (0.6- 
1), (0.7-1.1), and (0.7-1.3) for Via Lactea II, GHALO and GHALO s respectively. Thus, 
from the perspective of Poisson limits, these results would suggest that those derived from 
Maxwellian halos are possibly excessively aggressive by almost a factor of two. However, 
many limits are placed using one of Yellin's techniques [69] , where not only the overall rate, 
but also the distribution of signal versus background is important. Here we find that the 
halo uncertainties can be at their largest. We show in Fig. 9 the spectra at DAMA and 
CRESST for a 100 GeV WIMP with <5 = 130, 150, 170 keV. We employ energy smearing at 
DAMA by assuming that the smearing reported for the one of 25 targets of DAMA/LIBRA 
[70] is characteristic of all of them. We assume a smearing of 1 keV at CRESST. 

One can see that the peak positions and properties can change by quite a large amount. 
At DAMA, the effect of this is principally to shift the peak. Such an effect is largely 
degenerate with the quenching uncertainty, which can reasonably range from q = 0.06 to 
q = 0.09. (For a lower quenching value, for instance, the peak will shift to lower energy 
in keVee.) The effects would be similar at XENON10, where the energy smearing will 
eliminate most interesting structures, and the shift in the location of the peak will be 
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Figure 9: Left: The differential recoil energy distribution (in dru=events/kg/day/keV) at DAMA 
for selected spheres from VL2, GHALO and GHALO s , normalized to unity in the DAMA signal 
range 2-6 keVee Right: The distribution at CRESST for the same spheres, normalized to unity from 
0-100 keV. The models are m x = 100 GeV and S = 130 (top), 150 (middle), 170 (bottom) keV. 

comparable to the uncertainties induced from L e g- . In contrast, the spectrum for CRESST 
can vary dramatically, with the peak location moving from below 30 keV to above 60 keV. 
The implication of this is that techniques such as optimum interval, maximum gap and 
Pmax axe likely all overly aggressive in that a peak in the spectrum might arise at a specific 
location in the real Milky Way halo, but that would not be reproduced at the appropriate 
position, for instance in a Maxwellian halo. In general, any technique that relies on knowing 
precisely the predicted spectrum is unreliable when considering these variations. 

Although we are still limited by our sample of simulations, we can study the effects on 
various limits by looking at the allowed parameter space and limits in a variety of halos. 
To do so, we largely follow [27] in calculating the allowed parameter space and limits, with 
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Figure 10: As in Fig. 9, but with m x = 700 GeV. 



a few important differences. As pointed out by [25], the uncertainty in quenching factor 
for iodine can be extremely important. We thus consider the range q = 0.06 to q = 0.09 
and take the 90% allowed region to be the union of 90% allowed regions, and similarly for 
the 99% region. We smear the signal as described above. 

For CDMS we employ the maximum gap method with the data set specified in [27]. 

For limits arising from XENON10, we use the data presented in the recent reanalysis 
of [56], and calculate the maximum gap limit (as in ZEPLIN-II and ZEPLIN-III) for both 
[71] and [55] values of L e g, taking at every point the more conservative of the two. For 
limits, we employ the maximum gap method. Since the data are smeared and there are 
great uncertainties in L c g, the detailed spectral information arising from structure is lost. 

For ZEPLIN-II, we take the data described in [46] , but (conservatively) take the energy 
values to be shifted by a factor of 0.24/0.19 for a maximum gap analysis, which employs 
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Figure 11: Allowed parameter space for DAMA at 90% (purple) and 99% (blue) with m x — 
70 GeV. For comparison the 90% and 99% regions for DAMA/LIBRA rates only (as opposed to 
DAMA/Nal + DAMA/LIBRA) are shown in red and green lines. Constraints are CDMS (solid), 
ZEPLIN-III (long dashed, thin), ZEPLIN-II (long-dashed, thick), CRESST (medium-dashed) and 
XENON10 (short-dashed). The regions are MB with v = 220 km/s and v csc = 550 km/s (top left), 
VL2 (top right), and a sample sphere from VL2 (bottom, left) and GHALO s (bottom, right). 

the value of L e g from [71] at higher energies. 5 

For ZEPLIN-III, we utilize the data within the published acceptance box up to 15 
keVee for the maximum gap analysis. We additionally employ data outside the blind 
acceptance box up to 30 keVee, the maximum to which efficiencies were provided in [47]. 
We use the delineated la region for the data in this high energy range since the 3a region 
is not specified. To convert from keVee to keVr, we use 

E r = (0.142£ ee + 0.005) exp(-0.305 E°J m ) (4.2) 

We could parameterize L e g as a function of energy, and use different shifts at different energies, but 
the limits are essentially the same, being dominated by the events at the highest energy. 
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Figure 12: As in 11, but with m x = 150 GeV. 



below 10 keVee [25, 42], and a constant ratio of 0.55 above that. (We use a larger value 
than [42] at high energies to account for the energy dependent value of L c g- suggested by 
[71], which tends to weaken limits slightly.) Note that we do not employ a background 
subtraction using the science data as in [72], nor do we do so for ZEPLIN-II as the same 
basic technique seems to have led to a significant overprediction of background for ZEPLIN- 
III. 

For all Xe based experiments, we take a fixed smearing of \/30 keVr, which is typical 
of the smearing of the experiments in the range of the peak signal. 

For the CRESST experiment, we consider the data from the commissioning run [48], 
as well as the additional events at ~ 22 keV, 33 keV and 88 keV presented for Run 
31 by [73], and take a smearing of 1 keV. As described above, Yellin-style techniques 
are ironically inappropriate for experiments with high energy resolution. Because of the 
significant variation in the possible spectrum when varying halo models with the high 
energy resolution, we use a binned Poisson analysis, rather than a maximum-gap technique. 
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Figure 13: As in 11, but with m x = 300 GeV. 



We divide up the signal region into a low region {Er < 35 keV), a middle region (35 keV < 
Er < 50 keV) and a high region (50 keV < Er < 100 keV). In this way, we do not 
overweight the events near the zero of the form factor (in the mid region) , but still include 
some broad spectral information. We require that the rate be lower than 95% Poisson 
upper limits for each bin, which is a 95% limit when the signal is exclusively in the low bin 
and ~ 90% when there is signal in both the low and high bins, and ~ 85% in the (rarer) 
instances when there is signal in all three bins. Although stronger limits can be set using 
for instance, optimum interval, because of the uncertainties in the halo distribution, this 
seems inappropriate. In other words, if some of the events at CRESST are real, then the 
distribution of events may well be telling us the properties of the halo. 

We show the results of these limits in Figs. 11-14. Even when accounting for smearing, 
and using a binned Poisson limit instead of optimum interval, it is remarkable how signifi- 
cantly the allowed parameter space can change from simulation to simulation, and between 
spheres in the same simulation. We see that CRESST remains the most constraining, as 
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Figure 14: As in 11, but with m x = 700 GeV. 



found previously, but we see important differences. First, agreement is generally improved 
in these spheres. This is due to the presence of structure at high velocities which increases 
the modulated fraction at DAMA. As a consequence, small regions of allowed parameter 
space exist at high mass, the size of which depends significantly on the particular halo. 
Larger regions exist at smaller mass (~ 150 GeV). Moreover, the allowed region of param- 
eter space can shift dramatically in a and 5, with pockets appearing at large 5 from high 
velocity particles. Should iDM be relevant for nature, the detailed nature of the halo is 
clearly significant in determining the properties of the direct detection signals. 

4.2 Light dark matter 

Light (~ GeV) dark matter [5] has been proposed as a means to reconcile DAMA with the 
other existing experiments [6, 74, 51]. The current iteration of LDM involves "channelling" 
[75, 50] whereby for some small fraction of the time the entirety of the nuclear recoil energy 
is converted to scintillation light. Since the observed energy is lower, it allows lighter 
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Figure 15: Spectra for a sample of spheres best fit to the DAMA data (shown) within the different 
simulations for LDM models with m x — 3 GeV (top, left), 7 GeV (top, right), and 13 GeV (bottom). 



WIMPs to scatter at the DAMA experiment. Such light WIMPs may be incapable of 
depositing energy above the threshold at XENON10 or CDMS-Si, for instance, allowing 
DAMA to evade these bounds. 

Such an explanation is not without difficulty, however. These light WIMPs have a 
rapidly falling event spectrum, due to the decreasing probability of channeling at higher 
energies, and the increasingly suppressed number of particles with higher kinetic energy. As 
a consequence the spectrum of LDM seems generally to fall too rapidly at light masses [52, 
24, 53]. At higher masses, the spectrum is acceptable, but is excluded by other experiments. 
We should note again here that this is under the assumption that the errors in the data are 
statistical only. Should there be, e.g., significant changes in the efficiency, it is possible that 
LDM could give a better fit. As noted earlier, the uncertainties in L e g- should weaken these 
limits, while the most recent analysis from XENON10 (with no events to the S2 threshold) 
would strengthen the limits. Thus, while the models do not seem to describe the data well, 
there are adequate uncertainties that this scenario remains an interesting possibility as an 
explanation of the DAMA signal. 

Since these models are also extremely sensitive to the presence of high velocity particles, 
one might wonder whether, just as in the iDM case, the properties of halos at high velocities 
might modify the spectrum and the relative signal strength at other experiments. 

We begin again by addressing the question of the effects on the spectrum. We show 
these for three different masses in Fig. 15. One can see that, while there are noticeable 
changes in the spectrum, its overall shape remains basically unchanged. The difference 
between LDM and iDM here is simple: in iDM, the "threshold" (i.e., minimum velocity) 
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Figure 16: The allowed parameter space for LDM models. Top, left: MB with vq = 220 km/s and 
v csc = 550 km/s. Top, right: VL2. Middle, left: GHALO. Middle, right GHALO s . Bottom two 
sample spheres taken from the simulations. 

scattering is actually at large energy, and to go to lower energies one requires higher ve- 
locities. Competing against this are the significant form factor effects. The competition 
between these two sensitive quantities leads to dramatic changes and features. For LDM, 
the threshold scattering is at zero, and all effects (channeling, form factor, velocity dis- 
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tribution) contribute to the falling (unmodulated) spectrum. Thus, it can fall somewhat 
more or less rapidly, depending on the number of particles at high velocity, but without 
very significant changes. 

To determine the spectrum, we follow the analysis of [52]. The most proper analy- 
sis would be to include the updated measurement of L e g (which will weaken the limits) 
and the reanalysis to remove multi-hit events from [56] (which strengthens the limits). 
However, incorporating the S2 threshold, resulting in contribution to an energy-dependent 
acceptance at low energies is complicated. Lacking detailed information on this, we will 
restrict ourselves to the previous XENON10 limits. Since our point is the relative change 
from halo model to halo model, this is reasonable, so long as we recognize that the limits 
should be taken with a grain of salt. 

We show in Fig. 16 the allowed parameter space for the LDM scenario for the canonical 
MB model, for the (averaged) simulations VL2, GHALO and GHALO-scaled, as well as 
two sample spheres from the halos. We see that the small spectral effects translate into 
small effects in the allowed parameter space. There can be some improvement, however, 
and although the exclusion curves technically cover the allowed regions, uncertainties in 
the low-energy scintillation of Xe, for instance, and the ambiguity of employing statistical 
errors without a covariance matrix suggest that a small region may yet be allowed near 
10 GeV. We should again emphasize that there is no sense in which the sampling we have 
done is statistically complete, and our failure to find an allowed region does not preclude 
the possibility that one exists, or might arise in another simulation. A proper reanalysis 
using the data from [56] is warranted. 

5. Discussion 

Direct detection of dark matter intimately ties up questions of astrophysics and particle 
physics. In recent years, the range of particle physics models has exploded, many with new 
interactions and properties. Some of these, in particular light dark matter and inelastic dark 
matter, are particularly sensitive to the high velocity tail of the WIMP velocity distribution. 
The same is true, even for standard WIMPs, for many directionally sensitive detection 
experiments (e.g. DRIFT, NEWAGE, DMTPC), which rely on high recoil energies in 
order to measure the direction of the scattering DM particle. It is especially at these high 
velocities that deviations from the standard isothermal Maxwell-Boltzmann assumption 
are expected to be most significant. Velocity substructure, arising from nearby subhalos 
and tidal streams, as well as from the host halo's velocity anisotropy, can have profound 
effects on the expected event rates, the recoil energy spectrum, and the typical direction 
of scattering DM particles. 

We have studied the effects of these deviations by employing data from N-body simula- 
tions directly, rather than with analytic fits that miss interesting structures. Our analysis 
is based on two of the highest resolution numerical simulations of Galactic dark matter 
structure, Via Lactea II and GHALO. We confirm previously reported global departures 
from the Maxwell-Boltzmann distribution, with a noticeable excess of particles on the high 
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velocity tail. In many of our local samples we find additional discrete features due to 
subhalos and tidal streams. 

These global and local departures from the standard Maxwellian model have a number 
of interesting consequences: 

• For directionally sensitive detection, we have found that the fraction of particles 
coming from the "hottest" direction in the sky increases with the velocity threshold. 
More importantly, the direction from which most DM particles are coming can be 
affected significantly. At high velocity thresholds (> 500 km/s in the Earth's frame) 
the typical direction is more often than not shifted by > 10° away from the direction 
of Earth's motion, and up to 80° in extreme cases. 

• For inelastic DM the effects are dramatic. The relative strength of CDMS limits 
versus DAMA can change by an order of magnitude, while the relative strength of 
CRESST limits can change by almost a factor of two. For CRESST most of this 
effect is due to an enhanced modulation, which can be a factor of two larger than for 
a Maxwellian halo. Such effects would persist at Xenon detectors as well, changing 
their sensitivity relative to DAMA by 0(1). 

• For light DM, we have seen that small changes in the spectrum are possible, although 
not so much as to allow very light (~ 1 GeV) WIMPs to fit the existing data, but 
only somewhat heavier (^ 10 GeV). 

• In general, dramatic variations in the spectrum can appear for iDM models, compli- 
cating attempts to employ Yellin-type analyses, which rely upon a detailed knowledge 
of the predicted spectrum of the model. Instead we advocate a binned Poisson anal- 
ysis, which is less sensitive to assumption about the shape of the signal spectrum. 

• A further important effect that we have found is the possible variation of the phase 
of modulation as a function of energy. This has been noted before in the context 
of streams [76]. Since most background modulations would be expected to have the 
same phase as a function of energy, such a change could be a strong piece of evidence 
that a modulated signal was arising from DM. 

It is important to keep in mind that our simulations are nowhere close to resolving the 
phase-space structure at the relevant sub-parsec scales, and we rely on extrapolation from 
a coarser sampling (1-1.5 kpc radius spheres). If velocity substructure does not persist on 
smaller scales, then our analysis may overestimate the likelihood of these effects. On the 
other hand, we are probably underestimating the significance of the effects, since at the 
moment we are diluting the substructure signal by the host halo particles, while it would 
completely dominate the background host halo, should the Earth be passing through one 
of these substructures. 

We have seen here how, should DM be conclusively discovered with direct detection 
experiments, it may begin to help unlock questions about the formation of the galaxy, 
precisely because of these dramatic sensitivities. With the fantastic improvements in direct 



- 28 - 



detection on the horizon it is more important than ever to increase our awareness of the 
significant impact that the detailed structure of our Galaxy's dark matter halo may have. 
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